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1.  INTRODUCTION 


In  recent  years,  considerable  effort  has  been  devoted  to  the  aerodynamics 
of  vertical-takeoff-and-landing  (VTOL)  aircraft  in  ground  effect.  In  this 
mode  of  operation,  which  is  associated  with  takeoff,  hover,  and  landing,  there 
is  a  strong  interaction  between  the  lift-jet  flow  and  the  surrounding  fuselage 
and  ground  surfaces.  The  resulting  interference  forces  affect  significantly 
the  aircraft  performance.  Because  of  the  increased  complexity  of  the  ground- 
effect  flowfield,  which  is  characterized  by  three-dimensional,  turbulent, 
interacting  free  jets  and  wall  jets,  these  forces  cannot  be  predicted  by 
analytical  procedures  developed  to  describe  the  external  aerodynamics  of 
conventional-takeoff-and-landing  aircraft  without  extensive  modification. 

Attempts  to  predict  these  flows  currently  rely  on  a  so-called  modular 
analysis,  such  as  the  techniques  described  in  References  1  and  2.  With  these 
approaches,  the  complete  flowfield  is  divided  into  components,  specifically 
the  free  jets,  the  wall  jets,  and  the  upwash  formed  through  the  collision  of 
the  opposing  wall  jets.  The  components  are  modeled  using  a  separate  system  of 
equations  for  each  zone.  The  solutions  are  then  iteratively  pieced  together 
to  provide  a  description  of  the  total  region  of  interest.  One  shortcoming  of 
this  approach  is  the  difficulty  frequently  encountered  in  accurately  patching 
together  the  component  solutions.  In  addition,  the  modular  analysis 
techniques  rely  heavily  on  experimental  data  which  are  frequently 
configuration  dependent. 

In  an  effort  to  describe  jet  impingement  flows  from  a  more  general  and 
fundamental  point  of  view,  MDRL  has  performed  work  under  ONR  support  on 
solution  of  the  Reynolds-averaged  Navier-Stokes  and  turbulence  model  equations 
applied  to  these  flowfields.  Initially,  calculations  were  made  for  an 
isolated  planar  lift  jet  discharging  from  a  simulated  fuselage  undersurface  in 
ground  proximity  for  both  incompressible  and  compressible  conditions.  This 
work,  which  is  reported  in  Reference  3,  provides  accurate  predictions  of  the 
fuselage  undersurface  and  ground-plane  pressure  distributions,  illustrating 
their  variation  with  Reynolds  number  and  jet  height  above  ground.  In 
addition,  the  experience  gained  in  the  development  of  the  flowfield  model  and 
the  numerical  solution  scheme  for  two-dimensional  flow  formed  the  basis  for 
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extending  the  analysis  to  the  considerably  more  complex  three-dimensional  jet 
impingement  configurations. 

The  geometry  of  interest  in  the  present  study  consists  of  two 
interacting,  initially  axisymmetric  jets  with  fountain  formation;  the  geometry 
is  shown  schematically  in  Figure  1  for  the  general  case  of  inclined 
impingement.  The  jets  exit  from  two  axisymmetric  nozzles,  each  with  diameter 
D,  whose  centerlines  are  separated  by  a  spacing  S.  The  distance  from  the 
nozzle  exit  plane  to  the  midpoint  between  the  intercepts  of  the  geometric  jet 
centerlines  with  the  ground  plane  is  denoted  by  Hc ,  and  the  jet  inclination 
angle  relative  to  the  ground  plane  is  denoted  by  a. 
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The  impingement  of  the  two  parallel  jets  on  the  ground  creates  two 
adjacent  wall  jets  which  spread  radially  and  collide  to  form  a  fountain  upwash 
flow  that  is  directed  away  from  the  ground.  The  flow  properties  of  the 
fountain  are  determined  by  the  geometric  parameters  of  the  configuration  (S, 
Hc,  and  a)  and  the  relative  strengths  of  the  two  primary  jets.  For  large 
values  of  nozzle  centerline  spacing  (typically  S/D  >  10),  the  upwash 
characteristics  can  be  predicted  from  the  wall-jet  properties,  which  in  turn 
can  be  correlated  with  the  nozzle  exit  flow  conditions.  However,  for  closer 
centerline  spacing,  this  may  not  be  the  case  since  the  fountain  flow  can  be 
affected  by  its  interaction  with  the  primary  jets.  This  situation  becomes 
more  pronounced  as  the  jet  centerline  spacing  S  and  height  above  ground  H  are 
decreased. 

In  VTOL  aircraft  design,  close  nozzle  spacing  is  of  interest  since  it  can 
reduce  weight  requirements.  Consequently,  dual-jet  impingement  with  fountain 
formation  for  S/D  <  10  is  a  relevant  research  problem.  An  approximate  method 
for  modeling  close-jet  phenomena  using  a  modular  approach  based  on  empirical 
data  and  global  conservation-law  constraints  has  been  reported  in  Reference  4 
for  the  range  2  <  S/D  <  6.  The  flow  models  for  the  various  zones  are  coupled 
through  their  dependence  upon  the  upstream  behavior  of  the  prior  region  and 
the  downstream  position  of  the  opposing  jet.  The  technique  provides  the  flow 
properties  of  engineering  interest,  but,  as  pointed  out  by  the  authors 
(Refernce  4),  further  refinement  requires  a  more  extensive  experimental  data 
base  and  a  better  resolution  of  the  influence  of  turbulence-induced  static 
pressure  variations. 

The  analytical  approach  adopted  in  the  current  study  for  the  dual-jet 
impingement  analysis  is  a  solution  of  the  time-averaged  Navier-Stokes 
equations,  with  a  complementary  turbulence  model,  for  steady,  incompressible 
flow.  The  conservation  equations  are  solved  in  terms  of  scalar  and  vector 
potential  functions  and  vorticlty  using  a  second-order-accurate  finite- 
difference  algorithm.  Flowfield  characteristics  gained  from  such  a  study  can 
be  used  to  evaluate  and  improve  engineering  models  of  the  fountain  flow 
(References  1,  2,  and  4)  that  are  the  basis  for  predicting  VTOL  aircraft 
performance  in  ground  effect. 
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This  report  describes  the  flowfield  model,  the  numerical  solution 
technique,  and  representative  results.  The  latter  include  flowfield 
properties  for  the  following  dual-jet  impingement  configurations,  which  are  of 
interest  in  the  current  study:  equal-strength  jets  with  normal  impingement, 
unequal-strength  jets  with  normal  impingement,  equal-strength  jets  with 
inclined  impingement,  and  unequal-strength  jets  with  inclined  impingement. 

The  direction  of  future  MDRL  work  in  the  calculation  of  three-dimensional  jet 
impingement  flowfields  under  ONR  contract  N00014-80-C-0454  is  discussed. 


2.  THE  FLOWFIELD  MODEL 


In  this  section,  the  governing  equations  are  presented  which  are  used  to 
describe  the  turbulent,  incompressible,  steady  flow  of  the  three-dimensional 
impinging  jet  configurations  of  interest.  The  complete  time-averaged  Navier- 
Stokes  equations  and  a  one-equation  turbulence  model  are  given  in  terms  of  the 
velocity  components  and  pressure,  and  then  they  are  rewritten  in  terms  of 
vorticity,  scalar,  and  vector  potential  functions.  Boundary  conditions  are 
defined  for  the  interacting  jets  with  fountain  formation. 

2. 1  The  Governing  Equations 

The  time-averaged  continuity  and  momentum  (Navier-Stokes)  equations  for 
steady  incompressible  flow  are  given  below  in  tensor  notation  for  a  space- 
fixed  reference  through  which  the  fluid  moves. 

Conservation  of  mass: 


Conservation  of  momentum  in  the  ith  direction: 


(2-1) 


3  (u1u1) 
3x. 


+  iE _ L-  (7  - 

+  8Xi  3Xj  ^  ij 


P  u. 


V  ■ 


(2-2) 


The  conventional  notation  is  used  in  which  p  is  the  density,  u^  the  velocity 
component  in  the  ifc^  direction,  p  the  static  pressure,  and  Tjj  the  shear 
stress  in  the  itl*  direction  on  a  surface  normal  to  the  jc^  direction. 

The  previous  equations  were  derived  using  Reynolds  decomposition  in  which 
a  general  flow  variable  $  is  expressed  as  the  sum  of  a  time-averaged 
component  $  and  a  fluctuating  component  $  -  <p  +  The  result  is  a 

system  of  time-averaged  conservation  equations  which  have  the  same  form  as  the 
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instantaneous  conservation  equations  with  the  exception  of  the  Reynolds  stress 


term  -p  u^'u^ '  . 

In  order  to  obtain  a  closed  system  of  equations,  the  Reynolds  stress  term 
must  be  related  to  the  time-averaged  quantities  that  describe  the  mean  flow. 

In  the  present  work,  this  is  accomplished  using  the  one-equation 
incompressible  turbulence  model  of  Glushko,  Reference  5.  With  this  approach, 
the  effective  shear-stress  term  in  Equation  (2-2)  is  represented  by 


P  Ui  Uj 


2 

3 


Vk> 


(2-3) 


where  is  the  effective  viscosity  (the  sum  of  the  molecular  and  turbulent 
components,  Peff  “  U  +  pt)  and  k  is  the  turbulent  kinetic  energy.  The  latter 
is  computed  from  the  following  transport  equation: 


3(u.k) 
3x  . 


3ui  2  3ul 

3XJ  "  3  6ljPk  3XJ 


3k 

ue(Xr)  — 


+  U 


32k 
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(2-4) 


In  the  Glushko  model,  the  turbulent  component  of  the  effective  viscosity  is 
given  by 


ut  -  lie(r),  (2-5) 

where  e(r)  is  computed  from  the  relation 

e (r)  =  H(r)ar.  (2-6) 

In  Equation  (2-6),  r  is  the  Reynolds  number  of  the  turbulence  based  on  the 
length  scale  of  the  turbulence  L, 
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-  /~k  L 

r  = - 


(2-7) 


a  is  a  constant  with  value  0.2;  and  H(r)  is  a  function  defined  by 


1  1.25  <  —  <  « 


L  ro  J 

where  rQ  «  110.  With  the  turbulent  kinetic  energy  computed  from  Equation  (2- 
4) ,  the  turbulent  component  of  the  effective  viscosity  can  be  evaluated  using 
Equations  (2-5)  through  (2-8).  The  effective  viscosity  in  turn  determines  the 
shear  stress.  Equation  (2-3),  appearing  in  the  time-averaged  momentum 
equation. 

Although  the  partial  differential  equations  that  describe  the  flowfield 
are  written  in  terms  of  u^,  they  are  not  solved  for  the  velocity  components. 
Rather,  the  equations  are  solved  in  terms  of  scalar  and  vector  potential 
functions  and  vorticity  using  an  approach  proposed  by  Aregbesola  and  Burley, 
Reference  6.  In  this  technique,  which  is  analogous  to  the  stream- 
function/vorticity  approach  frequently  adopted  in  two-dimensional  problems, 
the  pressure  is  eliminated  from  the  momentum  equations  in  the  formulation  of  a 
vorticity  transport  equation.  The  pressure  field  is  then  computed  from  a 
Poisson  equation,  also  derived  from  the  component  momentum  equations. 

Consider  first  the  equation  for  the  conservation  of  mass.  To  identically 
satisfy  the  continuity  condition,  Equation  (2-1),  a  scalar  potential  *  and  a 
vector  potential  A^  are  introduced  which  are  related  to  the  velocity  field  by 
the  following  expression: 


3Aj 
'ijk  3x 


(2-9) 
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In  Equation  (2-9),  is  the  alternating  tensor.  The  scalar  potential  Is 

Introduced  to  facilitate  setting  the  boundary  conditions  on  the  velocity 
components.  Hirasaki  and  Heliums,  Reference  7,  have  shown  that  a  unique 
vector  potential  A^»  within  the  gradient  of  an  arbitrary  harmonic  function, 
exists  such  that 


—  -  0.  (2-10) 

3xi 


By  taking  the  divergence  of  Equation  (2-9)  and  introducing  the  conditions 
imposed  by  Equation  (2-1)  (the  conservation  of  mass)  and  Equation  (2-10),  a 
partial  differential  equation  for  the  scalar  potential  follows, 


324> 

SXjBx^ 


0. 


(2-11) 


A  differential  equation  for  the  vector  potential  is  derived  by  taking  the  curl 
of  Equation  (2-9)  and  using  Equation  (2-10),  which  gives 


^i 

3xj3xj 


-  n 


i’ 


(2-12) 


where  the  vorticity  is  defined  by 


3uk 

Qi  “  Eijk  3^* 


(2-13) 


Consider  next  the  equation  for  the  conservation  of  momentum.  Taking  the  curl 
of  Equation  (2-2)  and  utilizing  Equation  (2-13)  eliminates  the  pressure  from 
the  momentum  equation  and  results  in  a  transport  equation  for  the  vorticity. 
The  pressure  is  then  computed  from  a  Poisson  equation,  which  is  derived  by 
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taking  the  divergence  of  Equation  (2-2). 

Carrying  out  the  steps  outlined  In  the  previous  paragraph.  Introducing  a 
characteristic  length  D  and  characteristic  velocity  VQ  to  normalize  the 
equations,  and  expanding  the  equations  In  the  cartesian  coordinates  (x,  y,  2) 
with  velocity  components  (u,  v,  w)  results  in  the  system  of  dimensionless 
equations  given  below. 

Pois&wi  equation  for  the  scalar  potential: 


lii  ,  3^4  3^4 

3x2  3y2  3z2 


0 


(2-14) 


Poisson  equation  for  the  x  component  of  vector  potential: 


32A  32A  32A 
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3x  3y  3z 


n 

X 


(2-15) 


Poisson  equation  for  the  y  component  of  vector  potential: 


32A  32A  32A 


(2-16) 


Poisson  equation  for  the  z  component  of  vector  potential: 


2  2  2 

3  A  3ZA  34A 

- *  +  - L  +  - L 

2  2  2 
3x  3y  3z 


Q 

z 


(2-17) 
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Transport  equation  for  the  x  component  of  vorticity: 


Transport  equation  for  the  y  component  of  vorticity: 


10 


Transport  equation  for  the  turbulent  kinetic  energy: 


11 


Poisson  equation  for  the  static  pressure: 


3 v  3v  ,  3w  3_v  3_u  3w  3v  3v  j)w  3w 
3 y  3y  3y  3z  3z  3x  3z  3y  3z  3z 


(2-22) 


The  only  dimensionless  parameter  which  enters  the  system  of  equations  is 
the  Reynolds  number.  Re  ■  V0D/v.  The  location  in  the  flowfleld  where  VQ  and  D 
are  chosen  will  be  defined  subsequently  for  the  particular  configurations  of 
Interest.  The  Poisson  equation  for  the  scalar  potential  is  independent  of 
Reynolds  number  and  need  be  solved  only  once  for  a  given  geometry  and  set  of 
boundary  conditions. 
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2.2  The  Boundary  Conditions 


In  the  system  of  equations  presented  in  Section  2.1,  there  are  nine 
elliptic  partial-differential  equations  which  must  be  solved  for  the  following 
nine  flow  properties:  $,  A^,  Ay,  Az,  Qx,  ny,  flz,  k,  and  p. 

The  velocity  components  normal  to  the  surface  S  of  the  computational 
domain  with  volume  R  are  assumed  to  be  known  and  are  used  to  evaluate  the 
normal  derivatives  of  the  scalar  potential  on  the  boundaries: 


1± 

3n 


niui» 


(2-23) 


where  n^  denotes  a  unit  vector  normal  to  the  boundary.  Equation  (2-23) 
provides  the  necessary  boundary  conditions  with  which  to  solve  Equation  (2-14) 
for  $  within  an  arbitrary  constant.  Since  only  the  gradients  of  the  scalar 
potential  are  required  to  compute  the  velocity  components,  the  value  of  the 
constant  is  immaterial. 

The  boundary  conditions  on  Che  vector  potential  are  chosen  such  that 
the  normal  component  is  determined  by  the  solenoidal  condition.  Equation 
(2-10),  and  the  tangential  components  are  taken  to  be  zero.  For  a  flat 
surface,  these  boundary  conditions  are 


9A 

3^-°’  At  "  °»  <2"24> 


where  t  denotes  the  tangential  component.  Since  the  Poisson  equations  for  the 
components  of  the  vector  potential  are  dependent  on  the  vorticity,  which 
appears  as  a  source  term  on  the  right  side  of  Equation  (2-12),  they  cannot  be 
solved  independently  as  can  the  Poisson  equation  for  scalar  potential. 

Rather,  the  vector  potential  equations  in  combination  with  the  remaining 
equations  must  be  solved  iteratively  for  A  on  the  boundaries  as  well  as  in  the 
Interior  region  of  the  computational  domain. 
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Similarly,  the  vorticity  on  the  boundaries  must  be  computed 
iteratively.  At  each  cycle  in  the  iteration,  the  velocity  components  at  the 
boundaries  are  updated  from  Equation  (2-9), 


3$  9Ak 

Ui  °  ~  3x±  Eijk  3Xj’ 


(2-9) 


using  the  fixed  value  of  $  and  the  current  value  of  A^.  In  turn,  with  u^ 
updated  on  the  boundaries,  the  boundary  values  of  the  vorticity  are  recomputed 
from  Equation  (2-13), 


o 

i  “  Gijk  3Xj’ 


(2-13) 


and  used  during  the  next  iterative  cycle  in  the  solution  of  the  vorticity 
transport  equations.  The  turbulent  kinetic  energy  variation  is  related  to  the 
magnitude  of  the  vorticity  on  the  boundaries  through  a  simple  mixing-length 
formulation,  except  on  a  solid  surface  where  k  *=  0,  and  therefore  is  also 
computed  iteratively  during  the  course  of  the  calculations. 

The  solution  of  the  Poisson  equation  for  the  static  pressure  can  be 
deferred  until  a  convergent  solution  is  obtained  for  Ax,  Ay,  A^  Slx,  fty,  ft2, 
and  k.  The  boundary  conditions  required  for  Equation  (2-22)  follow  from  the 
pressure  gradients  given  by  the  component  momentum  equations. 

0£  _  2_ 

3x  Re 

+ 

(2-25) 
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3p.  _  _2 

3y  Re 


3e 

3 


(r)  / 3v  3u\  ?  3e (r)  3v  3e  (r)  /  3v  3w\ 
x  y3x  3yJ  3y  3y  3z  l  3  z  3y  I 


+  [1  +  e(r)] 


(3  2v  3  2v  3  2y\  4  3k  .  /  3v  3v  3v\ 

3x2  3y2  3zV  "33y  V3X  "3y  W  3V 


(2-26) 


3p  2 
3z  Re 


3e 


+  [1  +  e(r) 


(r)  /sw  3u\  3e  (r)  / £w  £v\  3e  (r)  3w 
x  1 3x  +  3z  )  +  3y  1 3y  3zi+Z  3z  3z 

..  /s2w  32w  3  2w\  43k  t/^w.  3w  3w\ 

)]  U2  3y2  a.2/  '  3  3z  V  3x  "9y  W  9z/ 


(2-27) 


The  right  sides  of  Equations  (2-25)  through  (2-27)  are  known  from  the  solution 
of  the  equations  for  the  scalar  potential,  vector  potential,  vorticity,  and 
turbulent  kinetic  energy. 

The  two  computational  domains  that  are  used  in  the  present  work  and  for 
which  the  normal  velocity  components  at  the  surface  must  be  specified  are 
shown  in  Figure  2,  where  all  dimensions  have  been  normalized  by  D.  For  the 
case  of  two  equal-  or  unequal-strength  jets  with  normal  impingement,  Figure 
2(a),  the  boundaries  z  -  0,  x  »  0,  and  z  «  S  are  jet  symmetry  planes.  The 
boundary  y  «  0  is  th«  ground  plane,  the  boundary  y  “  h  is  the  inflow  plane 
through  which  the  primary  jets  enter  the  computational  domain,  and  the 
boundary  x  *  w  is  the  outflow  plane  through  which  the  wall  jet  exits  the 
computational  domain.  For  the  case  of  two  equal-  or  unequal-strength  jets 
with  inclined  impingement.  Figure  2(b),  there  is  only  one  symmetry  plane, 
x  ■  0.  The  boundary  y  *  0  is  the  ground  plane,  the  boundary  y  -  h  is  the 
inflow  plane,  and  the  boundaries  z  -  0,  x  -  w,  and  z  -  L  are  outflow  planes. 
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The  manner  in  which  the  required  velocity  profiles  through  the  primary 
jets,  the  wall  jet,  and  the  fountain  are  computed  at  the  boundaries  of  the 
computational  domain  is  presented  in  Appendix  A  for  the  most  general  case  of 
interest,  unequal-strength  jets  with  inclined  impingement.  For  dual  jets  with 
normal  impingement,  the  simplification  that  can  be  introduced  in  the  analysis 
is  noted. 

In  the  specification  of  the  boundary  conditions,  it  is  necessary  that  the 
imposed  inflow  and  outflow  velocity  profiles  satisfy  the  global  conservation 
of  mass  for  the  computational  domain.  The  Poisson  equation  for  the  scalar 
potential.  Equation  (2-11),  has  a  solution  if  and  only  if 

f^hr1dV-  /H-s'  -°-  (2-28> 

R  s  s  s 

If  the  net  mass  flux  u^n^dS'  is  not  extremely  small  (typically  on  the 
order  of  lO-*’),  a  convergent  solution  of  Equation  (2-11)  cannot  be  obtained. 

In  the  present  work,  the  normal  velocity  components  required  in  the  analysis 
are  specified  such  that  the  trapezoidal  approximation  to  Equation  (2-28)  is 
satisfied  within  the  required  tolerance. 

The  finite-difference  technique  that  is  used  to  solve  the  coupled  system 
of  nonlinear  elliptic  equations  is  described  in  the  next  section. 
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3.  THE  NUMERICAL  SOLUTION  SCHEME 


This  section  is  a  discussion  of  the  finite-difference  techniques  used  to 
solve  the  partial  differential  equations  which  describe  the  flowfield.  The 
Poisson-type  equations  are  discretized  using  the  standard  central-difference 
algorithm,  and  the  transport-type  equations  are  discretized  using  an 
augmented-central-dif f erence  algorithm.  The  manner  in  which  the  coupled 
system  of  equations  is  iteratively  solved  and  the  convergence  characteristics 
of  the  solution  are  described. 

3. 1  The  Discretization  of  the  Poisson-Type  Equations 

The  Poisson  equations  for  the  scalar  potential,  the  three  components  of 
the  vector  potential,  and  the  pressure  can  all  be  written  in  the  following 
form  in  terms  of  an  arbitrary  flow  variable  ip  and  source  term  a: 


3x  3y  3z 


(3-1) 


Writing  Equation  (3-1)  for  the  scalar  potential,  $  is  replaced  with  $ ,  and 
a  Is  0,  which  reduces  the  equation  to  Laplace  form-  The  boundary  conditions 
are  fully  Neumann, 


3* 

3n 


niV 


(3-2) 


Writing  Equation  (3-2)  for  the  components  of  the  vector  potential,  <j>  is 
replaced  with  Ax,  Ay,  and  Az  for  a  given  by  -  f)x,  -  £1  ,  and  -  ft2, 
respectively.  The  boundary  conditions  are  mixed  Neumann  and  Dirichlet, 


3A 

=  0;  A±  =  0,  i  *  n.  (3-3) 
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Writing  Equation  (3-1)  for  the  static  pressure,  4>  is  replaced  with  p  and  o  is 
given  by  the  right  side  of  Equation  (2-22),  which  is  denoted  Op  for 
simplicity.  The  boundary  conditions  are  Neumann, 


1&  _ 

3n 


(3-4) 


where  the  pressure  gradient  normal  to  the  surface  is  given  by  either  Equation 
(2-25),  (2-26),  or  (2-27).  However,  in  the  solution  of  Equation  (3-3)  for  p, 
the  pressure  is  integrated  along  the  plane  x  =  0  of  the  computational  domain 
in  order  to  set  the  level,  which  replaces  one  of  the  Neumann  conditions  by  a 
Dirichlet  condition. 

To  write  Equation  (3-1)  in  finite-difference  form,  the  three-dimensional 
nodal  network  shown  in  Figure  3  is  introduced.  A  uniform  grid  spacing  h  is 
used  in  each  of  the  coordinate  directions,  so  that  x  =  h(i  -  1),  y  *=  h(j  -  1), 
and  z  *»  h(k  -  1),  where  i,  j,  and  k  are  the  respective  node  indices.  The 
conventional  central-difference  algorithm  is  introduced  for  which  the  second 
derivatives  centered  at  an  interior  point  in  the  nodal  network  are  given  by 


i.j.k 


i.J.k 


i.j.k 


h2 

*i.H-l.k  ~  2<^i,  1  ,k  +  *i.1-l,k 
h2 

*1.1. k+1  "  2*i.1,k  +  *i,1,k-l 

t_2 


(3-5) 


(3-6) 


(3-7) 


When  Equations  (3-5)  through  (3-7)  are  substituted  into  Equation  (3-1),  the 
discretized  form  of  the  Poisson  equation  is 
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*i+l,j,k  +  *i-l,J,k  +  ^i , j+1 ,k  +  +  *itJtk+l 


^i,j,k-l  ~  6*i,j,k  =  h2°i,j,k 


(3-8) 


Equation  (3-8)  applies  directly  to  the  interior  of  the  computational  domain: 

1  <  i  <  I,  1  <  j  <  J,  and  1  <  k  <  K,  where  I,  J,  and  K  denote  the  maximum 
values  of  the  respective  node  indices. 

The  manner  in  which  the  boundary  conditions  are  imposed  on  Equation  (3-8) 
depends  on  whether  they  are  Dirichlet  or  Neumann.  In  the  former  case,  if  a 
value  of  <j>  on  the  left  side  of  Equation  (3-8)  is  known,  it  is  transferred  to 
the  right  side  of  the  equation  and  is  eliminated  as  an  element  of  the  solution 
vector.  In  the  case  of  a  Neumann  boundary  condition.  Equation  (3-8)  is 
applied  at  the  boundary  to  solve  for  <j>  at  that  location.  The  value  of  <p  on 
the  line  outside  the  computational  domain  is  computed  by  treating  this  line  as 
an  image  of  the  line  adjacent  to  the  boundary  on  the  inside  of  the 
computational  domain.  The  image-line  value  of  $  is  related  to  the  interior¬ 
line  value  through  the  central-difference  discretization  of  the  Neumann 


I  =  maximum  \  index 


Figure  3.  The  three-dimensional  finite-difference  stencil. 
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boundary  condition. 

at  the  boundary  x  « 
given  by 


Consider 
h(I-l) . 


is  known 


,  for  example,  the  case  where 
The  finite-difference  form  of  this  derivative  is 


!± 

3  x 


I.J.k 


h+iihl  ~  *1-1, J.k 

2h  » 


from  which  it  follows  that 


(3-9) 


I+l,j,k  *1-1, j,k 


2h 

3x 


I.J.k 


(3-10) 


Equation  (3-10)  is  then  substituted  into  Equation  (3-8)  written  for  i  =  I. 
Neumann  conditions  on  the  other  boundaries  are  treated  in  an  analogous 
fashion. 

The  entire  system  of  algebraic  finite-difference  equations  can  be  written 
in  the  following  matrix  form: 

M  *  b  i  (3-11) 


where 


A  -  [Aj,  A2,  ...,  Aj 


and  A^  is  a  block  tridiagonal  matrix  with  diagonal  elements 


(3-12) 

(3-13) 


21 


The  superdiagonal  and  subdiagonal  elements  are  all  (6 ^ ,  except  the  first 
superdiagonal  element  and  the  last  subdiagonal  element,  which  are  (26^j). 

In  Equation  (3-12),  represents  a  coefficient  matrix  for  a  fixed  k 
plane.  The  terms  and  j,k-l  have  been  transferred  to  the  right 

side  of  Equation  (3-8) .  Gauss-Seidel  line  relaxation  is  used  to  solve 
Equation  (3-11),  which  requires  an  additional  transfer  of  terms  to  the  right 
side  of  Equation  (3-8).  For  example,  when  j  is  allowed  to  vary  in  a  fixed  k 
plane,  the  terms  and  j  R  are  included  in  the  right-side 

vector.  Specifically,  the  altered  form  of  b  becomes 


ai*J.k  ”  *1-1, j,k  “  *i+l,j,k  ~  ^i ,  j ,{ 


k  .  .  ,  +  2h«  -  2h«,  4  |± 

i.J.k-l  l.J  3y  1>j>k  J.J  3y 


)J 


(3-15) 


Elements  from  the  superdiagonal  and  subdiagonal  of  At,  are  included  in  b 

K  ik 

Further  details  about  the  solution  of  the  Poisson-type  equations  in 
conjunction  with  the  solution  of  the  transport-type  equations  are  given  in 
Section  3.3. 


3*2  The  Discretization  of  the  Transport-T-v 


qua t ions 


The  transport  equations  for  the  three  components  of  the  vorticity  and  the 
turbulent  kinetic  energy  can  be  written  in  the  following  form,  where  again  $ 
is  an  arbitrary  flow  variable  and  Oj  and  02  are  source  terms: 


W 


+  i!±  +  +  Re  B  |i  + 

3y2  3  z  y  3x 


3y 


+  Re6  ff  “  Reoi  +  °2  * 


(3-16) 


Writing  Equation  (3-16)  for  the  components  of  the  vorticity,  <f>  is  replaced 
with  ftx,  fty,  and  Sl2.  The  definitions  of  the  coefficients  a,  6,  y»  and  5  and 
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the  source  terms  Jj  and  0 2  follow  directly  from  an  inspection  of  Equations  (2- 
18)  to  (2-20).  The  boundary  conditions  are  fully  Dirichlet, 


(3-17) 


Writing  Equation  (3-16)  for  the  turbulent  kinetic  energy,  $  is  replaced  with 
k;  a,  B,  y ,  6,  oj,  and  a  2  follow  from  Equation  (2-21).  The  boundary 
conditions  are  mixed  Neumann  and  Dirichlet, 


3^  ‘  °*  k  *  °k 


(3-18) 


Equation  (3-16)  cannot  be  solved  for  an  arbitrary  Reynolds  number  using 
the  conventional  central-difference  approximations  to  the  derivatives.  The 
reason  for  this  difficulty  is  that  the  transport-type  equations  contain  the 
Reynolds  number  as  a  coefficient  of  each  convective  term.  Consequently,  with 
the  standard  central-difference  alogrithm,  the  discretized  system  of  equations 
is  diagonally  dominant  for  only  a  limited  range  in  Reynolds  number.  Diagonal 
dominance  is  necessary  to  obtain  convergence  in  the  iterative  solution  of  the 
discretized  system  of  equations. 

To  ensure  convergence  for  all  Reynolds  numbers,  the  second-order-accurate 
augmented-central-dif ference  (ACD)  algorithm  developed  by  Hoffman,  Reference 
8,  was  extended  to  the  three-dimensional  transport  equation  represented  by 
Equation  (3-16).  The  basis  of  this  method  can  be  illustrated  by  considering 
the  derivative  3$/3x  of  the  latter  equation.  At  an  interior  point  i,j,k  in 
the  nodal  network,  this  term  is  evaluated  in  terms  of  the  adjacent  points  in 
the  x  direction  with  the  following  truncated  Taylor-ser ies  representation  and 
standard  central-difference  approximation  to  the  first  derivative: 


3 

3x 


*i+1.1.k  -  *1-1.1. k  h^  3^4 

.  .  .  “  2h  "  6  3 

1 9  J  t  k  3  x 


i.J.fc 


j£  ll± 

"  51  3x5 


l.j.k 


(3-19) 
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3  3 

In  the  ACD  scheme,  the  derivative  3  $/3x  is  retained  and  expressed  in  terms 
of  lower-order  derivatives  by  differentiating  Equation  (3-16)  with  respect  to 
x.  The  result  is 


3<fr 

3x 


i.j.k 


Wl.k  -  ♦i-l.J.k  +  hf  Re  £ 

2h  6  K  " 


“  3x2 


i.J.k 


Re  S 


li.J.k 


+  0 


(Reh2) 


(3-20) 


where 


,  /_L!£l£  +  i!x!£  +  !Y 

lifj.k  \a  3x  3x  a  3x  3y  a  3x3y 


32£ 


1  li  i£  +  I  6  s2<fr  _  I  Hi) 

a  3x  3z  a  3x3z  a  3 xj  ^  j  r 


(3-21) 


The  derivatives  3<}i/3y  and  3<ji/3z  in  Equation  (3-16)  are  represented  in  an 
analogous  fashion  with  the  ACD  algorithm.  The  following  is  the  discretized 
form  of  the  transport-type  equation: 


Cli,j,k  ♦i+l.j.k  +  C2i,j,k  ♦i-l.J.k  +  C3i ,  j  ,k  *i,j+l,k 


C7i,j,k  ^i.j.k  "  Di,j,k 


(3-22) 


The  definitions  of  the  coefficients  are  given  below. 
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(3-23) 
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f  i.j.t 
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„  _1  /„  Re0_h  ,  Re2h2  $2>1 

*j*k  h2\  ~T~~ 


2  6  a 


'3i 


'41 


"?(• 


,  Re^h  ,  Re2h2  xl 

2  6 


t.. 


'5i,j,k  h2 


7  f 


2  6  a, 


^  Re5h  ,  Re2h2  52' 
2  +  6  ~a 


'6l,j,k  “  h2 


f' 


Re6h  ,  Re2h2  52' 

2  +  6  J, 


i.j.k 


i.j.k 


i.j.k 


7i.j,k  -  h2 


2  ^3«+^i!  +  MVli+ReV 

6  n  6  a  6 


i.j.k  +°2i,j,k  6 
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ReV 


*3 

+  YS3  +  6S5\ 
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i.j.k 


(3-24) 


(3-25) 


(3-26) 


(3-27) 


(3-28) 


(3-29) 


(3-30) 


Sli,j,k  13  defined  by  Equation  (3-21);  S31>Jfk  and  S51>j>k  are  defined  by  the 
following: 
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./iMli  +  ifl  324>  .  ±  9y  ii  ,  1  35  3<4 
.j.k  ya  3y  3x  a  3x3y  a  3y  3y  a  3y  3z 


+  -i  5  rV- 
a  3y3z 


.1^ 

a  3yy 


i.j.k 


(3-31) 
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(3-32) 


S 


5i,j  ,k 


I  UL  !±  +  1  6  +  1  ii 

a  3  z  3  x  a  3x3z  a3z3y 


1_  32»  (  1  36  3d>  1  3ql 
a  ^  3y3 z  a3z3z  a  3z 


i.j.k 


The  boundary  conditions  imposed  on  the  transport-type  equation,  either 
Dirichlet  or  Neumann,  are  discretized  in  the  manner  described  in  the  previous 
section. 

The  algebraic  finite-difference  equations  can  be  written  in  the  following 
matrix  form: 

C  ♦  *  d  .  (3-33) 

Unlike  the  discretized  Poisson  equation  for  which  the  coefficient  matrix 
contains  constant  elements  resulting  from  the  discretized  Laplacian,  the 
coefficient  matrix  for  the  transport-type  equation  contains  elements  which  are 
functions  of  the  flow  properties  and  whose  values  must  be  updated  during  the 
course  of  the  iteration.  Equation  (3-33)  is  solved  using  Gauss-Seidel  line 
relaxation . 

Additional  description  of  the  solution  of  the  Poisson-type  and  transport- 
type  equations  is  given  in  the  following  section. 

3.3  The  Solution  Algorithm  for  the  Coupled  System  of  Equations 

The  sequence  in  which  the  solution  of  the  Poisson-  and  transport-type 
equations  is  carried  out  is  illustrated  in  Figure  A.  Initially  the  Laplace 
equation  is  solved,  and  the  values  of  $  are  stored  for  the  evaluation  of  the 
vector  potential,  the  vorticity,  and  the  turbulent  kinetic  energy.  The 
calculation  of  the  latter  is  carried  out  in  the  order  Ag,  Q z ,  Ay ,  fty,  Qx, 
and  k.  An  inner  iteration  loop  is  established  in  which  six  of  the  unknowns 
are  held  constant  and  the  seventh  is  calculated.  An  outer  iteration  loop  is 
completed  when  the  inner-loop  calculation  of  all  seven  variables  has  been 
carried  out.  Following  a  converged  solution  for  A^ ,  and  k,  the  Poisson 
equation  for  pressure  is  solved.  The  pattern  of  successive  line  relaxation 
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which  is  followed  for  all  the  variables  is  shown  in  Table  1.  For  example, 
with  regard  to  $  ,  for  a  fixed  x-y  plane,  the  variable  z  is  set  (<p ^  j  and 
4*1  ,j ,k-l  are  transferred  to  the  right  side  of  the  equation  being  solved). 


V  Computer  code  1 


>  Computer  code  2 


>  Computer  code  3 


figure  4.  Computing  sequence  used  in  the  solution  of  the  Ihmtield  equations. 


TABLE  1.  THE  PATTERN  OF  SUCCESSIVE  LINE  RELAXATION 
FOR  SOLUTION  OF  THE  FLOW  FIELD  EQUATIONS. 
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Within  the  x-y  plane,  the  variable  x  is  then  set,  which  results  in  a  line  of  <j> 
values  with  y  as  the  variable.  This  system  of  equations  is  solved,  and  the 
resulting  values  are  updated  using  successive  line  relaxation  as  x  is 
varied.  For  all  nine  flow  variables,  the  final  convergence  criterion  (the 
difference  in  magnitude  of  the  unknown  between  successive  iterations)  is  taken 
to  be  10-4. 

The  flowfield  calculations  in  the  present  work  were  done  using  the 
computers  of  the  Systems  Technology  Program  (STP)  -  System  Simulation  Center 
(SSC)  operated  by  the  McDonnell  Douglas  Astronautics  Company,  Huntington 
Beach,  for  the  Department  of  the  Army.  Jobs  were  submitted  from  a  Data  100 
terminal  in  MDRL  via  telephone  lines  to  the  STP  CDC  6400,  which  acts  as  an 
input/output  station  and  host  computer  for  a  CDC  7600.  The  jobs  were  batched 
on  the  7600,  which  is  a  65k  SCM  (small  core  memory)  -  256k  LCM  (large  core 
memory)  configuration.  Printed  output  was  generated  at  MDRL  via  the  Data  100 
printer,  and  tape  output  was  mailed  to  St.  Louis  for  use  with  the  MDRL 
computer  graphics  routines  on  the  McDonnell  Douglas  Automation  Company 
CYBER  175  computer. 

Typical  computing  tines  on  the  7600  are  1.5  min  each  for  the  solution  of 
the  Laplace  equation  for  the  scalar  potential  and  the  Poisson  equation  for 
static  pressure  (approximately  1000  iterations  each)  and  30-45  min  for  the 
solution  of  the  Poisson  equations  for  vector  potential  and  the  transport 
equations  for  vorticity  and  turbulent  kinetic  energy  (approximately  120  outer 
iteration  loops).  Representative  flowfield  properties  are  presented  in  the 
next  section. 
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k. 


THE  COMPUTED  FLOWFIELDS 


Under  the  present  contract,  turbulent  flowfields  Induced  by  Interacting 
impinging  jets  in  ground  effect  were  calculated  for  the  eight  configurations 
defined  in  Table  2.  To  illustrate  details  of  the  flow  property  distributions, 
the  computed  variables  (scalar  and  vector  potential  functions,  vorticity, 
velocity,  turbulent  kinetic  energy,  and  static  pressure)  were  displayed  in  the 
form  of  contour  plots  on  selected  planes  passing  through  the  computational 
domain  in  the  three  coordinate  directions.  To  provide  visualization  of  the 
overall  flow  pattern,  three-dimensional  pathlines  were  constructed  from  the 
calculated  velocity  components.  In  this  section,  representative  contour  and 
pathline  plots  are  presented,  and  comparisons  between  computed  and  measured 
flow  properties  are  made.  A  technique  for  generalizing  the  boundary 
conditions  in  the  solution  scheme  is  described. 

4 . 1  Parallel  Jets  with  Normal  Impingement 

For  the  case  of  normal  jet  impingement,  consideration  was  given  to 
both  equal-  and  unequal-strength  jets.  In  the  former  case,  flowfield 
solutions  were  carried  out  in  which  the  jet  height  above  ground  was  kept 
constant  (Hc  *  4)  and  the  centerline  spacing  was  varied  from  5  to  12 
(S  »  5,  6,  9,  and  12).  In  Figures  5  and  6,  contours  of  $ ,  Ax,  Ay, 

Az>  ,  u,  v,  w,  and  k  are  shown  on  the  midplane  of  the  computational 

domain,  y  =  h/2,  for  S  =  5.  The  clearest  interpretation  of  the  flow  pattern 

I  \BI  t  1.  IMPIM.IV.  JH  (  OSHM  RM IONS  I  KK.VI  H>  IN  IHt 
HIIVU1U  I)  <  OMPt  I  AIIONs. 
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Ail  length*  are  normalized  by  the  jet  nozzle  diameter,  and  the  Reynold*  number 
is  bated  on  properties  evaluated  at  the  exit  plane  of  nozzle  2 
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and  the  entrainment  into  the  fountain.  The  analogous  contour  plots  for  the 
configurations  with  S  =  6,  9,  and  12  are  qualitatively  the  same  as  those  shown 
in  Figures  5  and  6. 


(c)  x -component  of  velocity,  u 


z 


(d)  y -component  of  velocity,  v 


(f >  Turbulent  kinetic  energy 


Figure  6.  Flow-field  properties  on  the  plane  >  =  h  2  for  equal-strength  jets  with  normal  impingement  (S  =  5.  h  =  2. 


Although  the  planar  contour  plots  provide  detailed  information  about  the 
flow  variables  throughout  the  computational  region,  they  do  not  provide  a 
simple  visualization  of  the  overall  flow  pattern.  The  physical  charac¬ 
teristics  of  the  flowfield  are  best  illustrated  through  three-dimensional 
particle-pathline  (streamline)  traces.  These  plots  are  constructed  in  the 
following  manner:  after  obtaining  a  convergent  solution  of  the  Navier-Stokes 
equations,  the  three  components  of  the  velocity  vector  are  stored  as  a 
function  of  position  in  the  computational  domain.  A  point  is  then  selected  in 
the  flowfield,  and  a  particle  is  moved  along  the  local  velocity  vector  at  that 
point  over  a  small  time  increment.  When  the  particle  arrives  at  a  new  posi¬ 
tion,  it  is  again  directed  along  the  local  velocity  vector.  This  procedure  is 
repeated  stepwise  in  time  until  the  particle  exits  the  computational  domain. 
The  CRT  graphics  routine  has  the  capability  of  rotating  the  computational 
domain  about  each  of  the  three  axes  so  that  the  pathlines  are  visible  from  any 
desired  angle. 

Figure  7  illustrates  the  streamline  plots  for  the  case  of  equal-strength 
jets  with  normal  impingement  for  S  =  6,  9,  and  12.  The  origins  of  the  traces 
are  taken  to  be  the  nodal  points  at  the  periphery  of  the  primary  jets  on  the 
top  surface  of  the  computational  domain,  as  illustrated  in  the  rotated  ^op 
views  of  the  three  configurations.  The  views  at  the  left  in  Figure  7  are 
obtained  by  rotating  the  computational  domain  45°  counterclockwise  about  the 
z  axis.  Shaded  areas  have  been  added  to  the  top  plane  of  the  computational 
domain  to  define  the  primary  jet  and  fountain  regions.  These  plots  clearly 
show  the  change  in  the  entrainment  flow  into  the  fountain  as  S  is  doubled  in 
magnitude . 

The  accuracy  of  the  computations  with  regard  to  the  prediction  of  wall- 
jet  properties  for  the  case  of  normal  impingement  is  illustrated  in  Figure 
8.  The  experimental  results  are  for  a  single  turbulent  jet  striking  the 
ground  with  a  =  90°,  Reference  9,  and  the  computational  results  are  for  dual 
jets  with  normal  impingement  and  S  =  12.  For  this  value  of  centerline 
spacing,  such  a  comparison  is  reasonable.  As  illustrated  in  Figure  8(a),  the 
flowfield  near  the  wall  is  characterized  by  an  impingement  region  in  which  the 
flow  changes  direction  and  a  wall-jet  region  in  which  the  shear  layer  develops 
over  the  ground  plane. 
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A  velocity  profile  computed  in  the  wall-jet  region  is  shown  in  Figure 

8(b)  in  which  w/wmax  is  plotted  as  a  function  of  y/yj/2»  where  y^/2 

value  of  y  for  which  w/w  „  =0.5.  Outside  the  inner  turbulent  boundary  layer 

J  max  y 

(y  >5),  Abramovich's  relation  for  the  velocity  profile,  Reference  9,  is  appli¬ 
cable. 


~~  -  (i  -  c3/2)2, 

w 

max 


(4-1) 


where  £  =  (y  -  6)/b,  with  b  denoting  the  distance  between  the  point  of  maximum 
velocity  and  the  edge  of  the  shear  layer.  Equation  (4-1)  and  the  profile 
obtained  from  the  Navier-Stokes  equations  coincide. 

Figure  8(c)  shows  the  measured  variation  (Reference  9)  of  the  maximum 
wall-jet  velocity  with  distance  from  the  stagnation  point  for  two  values  of 
jet  height  above  ground,  Hc  «  4  and  7.  The  local  maximum  velocity  wmax(rw) 
for  a  given  distance  from  the  stagnation  point  is  normalized  by  the  largest 
value  of  wMX  for  the  entire  range  of  rw,  a  value  denoted  by  wpeak.  In  the 
vicinity  of  the  stagnation  point  (rw  less  than  approximately  1),  viscous 
effects  are  negligible  and  wraax(rw) /wpeak  varies  essentially  linearly  with 
rw.  As  rw  increases  beyond  unity,  viscous  effects  increase  and  the  maximum 
wall-jet  velocity  decays.  In  the  absence  of  dissipation,  the  curves  would 
continue  to  increase  to  an  asymptotic  value,  as  shown  by  the  inviscid  flow 
curves  which  were  obtained  from  Bernoulli's  equation.  The  variation  of 
wmax  ^rw^wpeak  8*ven  by  the  Navier-Stokes  equations  shows  a  trend  identical 
to  that  observed  experimentally. 

The  computed  pressure  distribution  on  the  ground  plane  line  x  •  0  is 
shown  in  Figure  9(a)  for  the  case  of  dual-jet  impingement  with  S  =  12  and 
Re  »  100.  Also  shown  on  the  plot  are  the  data  of  Jenkins  and  Hill,  Reference 
10,  which  were  measured  for  the  same  configuration  but  with  a  Reynolds  number 
on  the  order  of  10^.  The  jets  stagnate  on  the  ground,  deflect  90°,  spread  as 
radial  wall  jets,  and  then  collide  at  the  midplane  between  the  jet  centerlines 
to  form  the  upwash  flow.  The  measured  pressure  profile  shows  a  rapid  decay 
from  the  stagnation  point  value,  an  essentially  constant  ambient  level  through 
the  wall  Jet,  and  a  rise  in  the  collision  zone  as  the  flow  decelerates  in  the 
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fa)  Ground  pressure  distribution  for  equal-strangth  jets  with  normal  impingement, 
x  *  0  (S  *  12,  h  *  2,  w  »  3,  Re  *  1001 


z 


(b)  Flow  direction  in  the  midplane  of  the  fountain  for  equal -strength  jets  with  normal  impingement 
(h-2,w«3,  Re  -  1001 


tFlow  direction  measured  by 
Jenkins  and  Hill  (Reference  10) 


S  -  6 


Jet 

centerline 


w 


Figur*  9.  Comparison  of  computed  and  measured  flow  properties  for  equal-strength  Jets  with  normal  Impingement. 


horizontal  direction  and  turns  in  the  vertical  direction.  The  computed 
pressure  rise  does  not  demonstrate  the  same  decay  rate  from  the  peak  pressure 
nor  the  rise  in  the  collision  zone.  This  discrepancy  is  attributed  primarily 
to  the  inability  of  the  relatively  coarse  grid  used  in  the  finite-difference 
calculations  to  adequately  represent  the  velocity  gradients  in  the  impingement 
and  upwash  regions. 

Also  shown  in  Figure  9  is  a  comparison  of  the  computed  and  measured  flow 
directions  in  the  midplane  of  the  fountain  for  equal-strength  jets  with  S  =  6 
and  9.  Most  zonal  analysis  schemes  which  are  used  to  model  the  upwash  for 
interacting  jets  are  based  on  the  assumption  that  the  radial  flow  pattern  of 
the  wall  jets  on  the  ground  plane  is  continued  into  the  upwash,  with  the 
virtual  origin  taken  to  be  a  distance  S/2  below  the  ground.  The  flow 
direction  was  determined  by  Jenkins  and  Hill,  Reference  10,  using  an  array  of 
17  flags  in  the  plane  z  =  S/2  which  were  mounted  on  a  ladder-shaped  wire 
support  permitting  each  flag  to  pivot  in  the  x  direction.  Photographs  taken 
to  show  the  flow  direction  and  used  to  construct  the  arrows  in  Figures  9(b) 
and  (c)  indicate  that  the  assumption  of  radial  flow  is  a  good  one  except  in 
regions  far  from  the  central  upwash,  a  condition  which  Jenkins  and  Hill 
attribute  to  entrainment  and  turbulent  mixing  of  ambient  air  into  the 
upwash.  The  comparison  of  the  experimental  flow  direction  with  the  computed 
streamlines  shows  that  the  latter  give  a  good  representation  of  the  flow 
direction  except  near  the  outflow  plane  x  =  w  where  the  streamlines  are  forced 
to  intersect  the  surface  in  a  direction  set  by  the  imposed  boundary 
conditions . 


Particle  pathline  plots  computed  for  unequal-strength  jets  with  normal 

impingement  are  shown  in  Figure  10.  In  this  case  the  jet  centerline  spacing 

was  held  constant,  and  the  momentum  ratio  at  the  exit  plane  of  the  two  jets 

was  varied  over  the  range  M^/M2  =  1«0,  0.5,  and  0.25.  These  momentum  ratios 

correspond  to  jet  centerline  velocity  ratios  of  V.  /V.  =  2.0,  0.707,  and 

l  ■*  2 

0.5.  With  the  primary  jet  and  fountain  regions  shaded  in  the  rotated  top  view 
of  the  computational  domain,  the  pathline  plots  show  the  shift  in  the  fountain 
toward  the  weaker  jet  and  the  strong  interaction  between  the  latter  and  the 
upwash . 


The  location  of  the  stagnation  line  resulting  from  the  collision  of  the 
opposing  wall  jets  is  indicated  by  the  contours  of  the  z-component  of  velocity 
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velocity  ratio,  there  is  a  displacement  of  approximately  one  diameter  in  the 
w  *=  0  line  from  the  midpoint  between  the  jet  centerlines. 


Figure  II.  Contours  of  the  {-component  of  velocity  on  the  plane  >  =  h  8  for  unequal-strength  jets  with  normal 
impingement  (S  =  9,  h  =  2,  w  =  J,  Re  =  1001. 


4.2  Parallel  Jets  with  Inclined  Impingement 


For  the  case  of  inclined  jet  impingement,  flowfield  calculations  were 
carried  out  for  both  equal-  and  unequal-strength  jets.  Contour  plots  of  the 
computed  flow  variables  for  equal-strength  jets  with  S  =  4.5,  L  *  10.5,  and 
a  »  85°  are  shown  in  Figures  12  and  13  for  the  midplane  y  =  h/2  of  the 
computational  domain.  Jenkins  and  Hill,  Reference  10,  do  not  report  a 
measured  ground-plane  pressure  distribution  for  this  configuration,  but  they 
do  provide  an  experimental  ground-plane  stagnation  line  which  is  compared  with 
the  computed  line  in  Figure  1A. 

A  calculation  of  unequal-strength  jets  with  inclined  impingement  was  made 

for  a  configuration  with  S  =  8,  L  =  14,  a  =  80°,  and  V._„  /V,„„  *=  0.707. 

jcej  jce2 

Figure  15  shows  contours  of  the  y-component  of  velocity  on  three  x-z  planes  in 
the  computational  domain  (y  =  0.33,  1.67,  and  2)  which  demonstrate  the 
computed  variation  in  the  primary  jet  velocity  profiles  as  the  ground  plane  is 
approached . 


4.3  Planar  Flowfields  Computed  with  Coordinate  Transformations 

A  shortcoming  in  the  present  flowfield  prediction  scheme  for  three- 
dimensional,  interacting,  impinging  jets  with  fountain  formation  is  the  rather 
involved  procedure  for  specifying  boundary  conditions.  This  approach, 
outlined  in  Appendix  A,  requires  empirical  data  and  numerous  engineering 
approximations  regarding  the  characteristics  of  the  flow.  A  better  technique 
is  to  remove  the  computational  boundaries  from  the  near  field  and  place  them 
in  the  far  field  where  simple  conditions,  such  as  no  property  gradients  in  the 
primary  flow  direction,  can  be  meaningfully  imposed.  This  approach  can  be 
adopted  through  the  use  of  coordinate  transformations  to  map  a  relatively 
large  physical  domain  into  a  smaller  computational  domain. 

An  exploratory  study  of  transforming  the  coordinates  in  the  Navier-Stokes 
equations  was  investigated  by  Dr.  R.  K.  Agarwal  (MDRL)  under  the  present 
contract.  An  Euler  transformation  was  used  to  map  the  physical-domain 
variables  (x,  y)  into  the  computational-domain  variables  (x,  y) , 
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<b)  x-component  of  vector  potential,  Ax 
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Figure  14.  (  ompamon  of  measured  and  computed  stagnation 
line  shape  fur  equal-strength  jels  with  inclined  impingement 
(S  =  4.5.  H(  =  2.74.  Or  =  85  .  Ke  =  100). 


The  constants  P  and  Q  were  chosen  to  obtain  the  desired  clustering  of  grid 
points  in  different  regions  of  the  physical  domain.  For  simplicity,  the 
coordinate  transformation  scheme  was  applied  to  the  two-dimensional,  time- 
averaged  Navier-Stokes  and  turbulence  model  equations  described  in  Reference  3 
rather  than  to  the  three-dimensional  equations.  The  Poisson  equation  for  the 
stream  function  and  the  transport  equations  for  vorticity,  turbulent  kinetic 
energy,  and  turbulent  dissipation  were  written  in  the  transformed 
coordinates  (x,  y)  and  solved  using  the  augmented-central-dif f erence  (ACD) 
algorithm. 

Computations  were  carried  out  for  the  following  t  iree  cases:  laminar 
entry  flow  in  a  channel,  impingement  on  a  flat  plate  of  a  planar  laminar  jet 
with  a  free  upper  boundary,  and  impingement  on  a  flat  plate  of  a  planar 
turbulent  jet  with  a  free  upper  boundary.  Flow  properties  computed  for  these 
geometries  are  shown  in  Figure  16.  For  the  laminar  entry  flow  (Re  =  75),  a 
fully  developed  condition  was  assumed  100  channel  widths  downstream  of  the 
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Figure  IS.  (  onlours  of  the  >-component  of  velocit>  for  unequal-strength  jets  with  inclined  impingement  (S  =  8,  I.  =  14.  a  =  80°, 
h  =  2,  w  =  3,  Re  =  100). 


entrance  plane,  and  Equation  (4-2)  was  applied  with  P  =  99/100  and  Q  =  0. 
Figure  16(a)  shows  the  variation  of  the  centerline  velocity  along  the  channel 
axis  and  the  excellent  agreement  with  the  computations  of  Wang  and  Longwell, 
Reference  11,  For  the  laminar  impinging  jet,  the  solution  was  carried  out 
with  the  nozzle  exit  plane  located  a  distance  Hc/D  =  40  above  the  ground  plane 
and  with  Re  =  10^.  In  the  transformation  given  by  Equation  (4-2), 

P  =  Q  =  39/40.  The  surface  vorticity,  which  is  shown  in  Figure  16(b),  assumes 
a  constant  value  in  approximately  three  jet  widths  from  the  stagnation 
point.  For  the  turbulent  impinging  jet  solution,  Hc/D  =  31,  Re  =  5.3  x  10  , 
and  P  =  Q  °  30/31.  Figure  16(c)  illustrates  the  ground  plane  vorticity. 

The  coordinate  transformation  work  with  the  Navier-Stokes  equations  for 
planar  flow  has  demonstrated  the  feasibility  of  extending  the  procedure  to  the 
three-dimensional  impinging  jet  configurations  with  fountain  formation. 
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(a)  Velocity  variation  along  the  centerline  for  a  laminar  duct  flow 


(b)  Vorticity  variation  along  the  plate  in  the  impingement  region  for  a  lamina/  jet 


(c)  Vorticity  variation  along  the  plate  in  the  impingement  region  for  a  turbulent  jet 


x/D  (Distance  along  the  plate  from  the  stagnation  point) 

GP030849  16 


Figure  16.  Computed  properties  for  t Ao-dimensional  duel  and  jet  impingement  flows  obtained  using 
coordinate  transformations. 
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5 .  SUMMARY 


Starting  with  the  complete  elliptic  conservation  equations  for 
incompressible,  steady,  three-dimensional  viscous  flow,  a  flowfield  model  has 
been  formulated  and  solved  numerically  for  two  parallel  jets  which  impinge  on 
a  ground  plane  and  interact  to  form  a  fountain  upwash.  The  time-averaged 
Navier-Stokes  equations  in  combination  with  the  Glushko  turbulence  model  are 
cast  in  terms  of  Poisson  equations  for  scalar  and  vector  potential  functions 
and  static  pressure  and  in  terms  of  transport  equations  for  vorticity  and 
turbulent  kinetic  energy-  This  system  was  solved  using  finite-difference 
procedures  for  eight  two-jet  conf igurations  with  various  centerline  spacings, 
angles  of  inclination,  and  nozzle  exit-plane  velocity  ratios.  In  this 
section,  the  characteristics  of  the  numerical  scheme  and  computed  flowfields 
are  summarized,  and  recommendations  are  made  for  generalizing  and  improving 
the  accuracy  of  the  prediction  method. 

5.1  Conclusions 

With  regard  to  the  numerical  solution  scheme,  no  significant  problems 
were  encountered  in  obtaining  convergent  solutions  of  the  coupled  system  of 
equations  for  the  conditions  defined  in  Table  2.  The  computed  flowfields 
provide  reasonable  prediction  of  the  wall-jet  characteristics,  the  upwash  flow 
direction,  and  the  stagnation  line  patterns  on  the  ground  plane.  However,  the 
ground-plane  pressure  profiles  do  not  predict  the  experimental  variation  in 
the  stagnation  zones  of  the  impinging  primary  jets  and  the  colliding  wall 
jets. 

The  limitations  of  the  present  prediction  method,  which  are  common  to 
most  Navier-Stokes  steady-state  finite-difference  schemes  applied  to  complex 
three-dimensional  flows,  are  the  following: 

(1)  The  use  of  empirically  determined  boundary  conditions  in  the  near 

field.  In  the  present  work,  the  rather  involved  procedure  outlined  in 
Appendix  A  for  setting  the  ne~--field  boundary  conditions  is 
necessary.  This  approach  requires  the  use  of  empirical  data  to 
establish  the  velocity  profiles  through  the  wall  jets  and  fountain  on 
the  boundary  planes  of  the  computational  domain.  A  oetter  approach  is 
to  take  the  computational  boundaries  in  the  far  field  where  simple 
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constraints  are  valid,  such  as  no  gradients  in  the  fluid  properties  in 
the  outflow  direction. 

(2)  The  use  of  coarse  grids  imposed  by  computer  storage  constraints.  In 
the  present  work,  the  same  grid  spacing  h  (the  distance  between  nodal 
points  in  the  finite-difference  mesh)  is  used  in  all  three  coordinate 
directions.  Consequently,  for  a  given  computational  domain,  if  the 
grid  spacing  is  decreased  by  a  factor  of  two,  for  example,  to  obtain  a 
greater  resolution  in  one  coordinate  direction,  h  must  be  divided  by 
two  in  the  remaining  coordinate  directions.  This  procedure 
essentially  doubles  the  number  of  grid  points  in  each  coordinate 
direction,  which  means  increasing  by  a  factor  of  eight  the  number  of 
storage  locations  required  for  each  triply  dimensioned  variable.  In 
the  present  code  there  are  20  arrays  with  dimensions  (I,  J,  K)  with 
the  result  that  halving  the  grid  spacing  increases  the  dimensioning  of 
these  variables  by  a  factor  of  160.  Consequently,  an  improved 
approach  is  to  permit  specification  of  the  node  spacing  independently 
in  each  of  the  three  coordinate  directions.  Moreover,  there  is  a  need 
to  nonuniformly  cluster  the  grid  points  within  a  given  coordinate 
direction  so  that  the  regions  of  the  flow  where  property  gradients  are 
the  largest  (such  as  the  wall-jet  regions  near  the  stagnation  zones) 
can  be  computed  with  greater  resolution.  This  approach  should  also 
permit  solutions  to  be  carried  out  at  a  Reynolds  number  several  orders 
of  magnitude  larger  than  that  considered  in  the  present  study. 

5.2  Recommendations 

The  two  limitations  described  previously  can  be  removed  through  the  use 
of  ci-ordinate  transformations  to  permit  the  application  of  the  boundary 
conditions  in  the  far  field  and  to  permit  clustering  of  the  grid  points  in 
regions  of  large  property  gradients. 

With  regard  to  the  specification  of  the  boundary  conditions,  the 
objective  of  the  transformation  is  to  map  a  large  (even  infinite)  physical 
domain  for  which  the  boundary  conditions  are  known  into  a  smaller 
computational  domain  which  can  be  adequately  represented  by  a  finite  number  of 
grid  points.  With  regard  to  improving  the  resolution  of  the  flowfield, 
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coefficients  can  be  introduced  into  the  coordinate  transformations  which 
permit  clustering  the  grid  points  in  regions  of  large  property  gradients.  As 
pointed  out  by  Oberkampf,  Reference  12,  by  transforming  to  a  computational 
plane  which  "stretches"  regions  of  large  property  gradients  in  the  physical 
plane  and  "compresses"  regions  of  small  property  gradients  in  the  physical 
plane,  it  is  possible  to  reduce  the  truncation  error  of  the  finite-difference 
solution.  Generally,  the  computations  can  be  carried  out  more  efficiently  and 
for  a  higher  Reynolds  number  than  can  be  done  using  a  relatively  coarse 
uniform  grid. 

Under  contract  to  ONR,  MDRL  will  use  coordinate  transformations  in 
solution  of  the  three-dimensional  time-averaged  Kavier-Stokes  and  turbulence 
model  equations.  In  addition,  the  augmented-central-dif ference  (ACC) 
discretization  scheme  used  in  the  work  presented  in  this  report  will  be 
replaced  by  an  upwind  difference  scheme  to  eliminate  the  storage  of  the  triply 
dimensioned  arrays  a,  (? ,  y,  and  6  [Equation  (3-16)]  required  in  the  ACD 
algorithm.  The  upwind  difference  formulation  will  permit  more  grid  points  to 
be  used  to  better  resolve  the  flowfield.  Calculations  will  be  made  for  equal- 
and  unequal-strength  jets  with  normal  and  inclined  impingement. 
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APPENDIX  A:  DEFINITION  OF  THE  VELOCITY  BOUNDARY  CONDITIONS 
FOR  THE  THREE-DIMENS IONAT.  JET  IMPINGEMENT  CONFIGURATIONS 


This  appendix  presents  details  of  the  calculation  of  the  velocity 
profiles  normal  to  the  surfaces  of  the  computational  domain  for  the  most 
general  case  of  interest  in  the  present  analysis,  unequal-strength  jets  with 
inclined  impingement.  The  simplifications  in  the  specification  of  boundary 
conditions  for  the  case  of  normal  impingement  of  two  jets  are  noted. 


Consider  the  configuration  shown  in  Figure  17.  Both  jets  have  a  nozzle 
exit  diameter  D,  and  Hc  denotes  the  distance  from  the  nozzle  exit  plane  to  the 
midpoint  of  the  distance  between  the  centerline  intercepts  with  the  ground 
plane.  The  angle  of  inclination  of  each  jet  is  ajj»  and  the  nozzle  centerline 
separation  distance  is  S.  The  computational  domain,  indicated  by  the  dashed 
lines,  has  dimensions  w,  h,  and  L  in  the  x,  y,  and  z  directions, 
respectively.  For  the  definition  of  the  velocity  boundary  conditions,  the 
origin  has  been  shifted  from  the  left  boundary  to  the  center  of  the 
comoutational  domain. 


First,,  consider  the  specification  of  the  velocity  profiles  in  the 
entering  jets  at  the  top  surface  of  the  computational  domain,  y  «=  h. 
Bradbury's  profile  based  on  experimental  results,  Reference  14,  is  used  in 
which  the  velocity  variation  V  through  the  Jet  normalized  by  the  centerline 
value  Vjc  is  given  by 


-  exp  [-0.6749n ,  (1  +  0.0269nJ)],  (A-l) 

Jc  J  3 


where  *  r^/6,  rj  is  the  radial  coordinate  from  the  jet  centerline,  and  6  is 
the  value  of  rj  at  which  V  ■  0.5  vjc*  A  plot  of  this  profile  is  shown  in 
Figure  18(a).  Vjc  and  5  are  obtained  from  the  empirical  curve  fits  shown  in 
Figure  18(b)  of  Vjc/Vjce  and  6/y"  as  functions  of  y"/D.  vjce  is  the  nozzle- 
exit-plane  jet  centerline  velocity,  and  y"  is  a  coordinate  perpendicular  to 
the  nozzle  exit  plane. 
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Figure  17.  Definition  of  the  lion  field  parameters  for  two  jets  of  equal  or  unequal  strength  *ith  inclined  impingement 
on  a  ground  plane. 


In  the  configuration  of  Figure  17,  the  distance  from  the  exit  plane  of 
nozzle  1  to  the  point  where  the  geometric  centerline  intercepts  the  ground 
plane  is  given  by 


H 


1 


H 

c 


S  ctn  a 


11 


2 


(A-2) 
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Similarly,  the  corresponding  distance  from  nozzle  2  is 


S  ctn  a 


H  + 
c 


’ll 


(A-3) 


52 


The  geometric  centerlines  of  the  two  nozzles  intercept  the  top  surface  of  the 
computational  domain  at  the  following  distances  from  the  nozzle  exit  planes: 


v"  =  H  - — 

yh ,  1  1  sin  a 


JI 


y"  =  H - — 

yh,2  2  sin  a 


JI 


(A-4) 


(A-5) 


Using  the  data  of  Figure  18(b),  the  values  of  Vjc^  j,  6^  Vjc^  2«  and 

5  are  computed.  The  required  velocity  components  at  the  top  surface  of  the 
n ,  z 

computational  domain  are  then  evaluated,  completing  the  definition  of  the 
inflow  velocity  profiles  through  the  primary  jets. 

Second,  consider  the  specification  of  the  velocity  profiles  in  the 
exiting  wall  jets  at  the  side  surfaces  of  the  computational  domain, 

z  =  -  y  ,  z  =  \  «  and  x  =  w.  It  is  first  necessary  to  define  the  origins  from 

which  the  wall  jets  emanate.  As  pointed  out  by  Jenkins  and  Hill  in  Reference 
10,  the  impingement  of  an  inclined  jet  on  a  ground  plane  produces  a  wall  jet 
in  which  the  flow  properties  depend  on  both  radial  distance  from  the  origin 
and  azimuthal  orientation.  If  <f>  denotes  the  angular  measurement  around  the 
jet  impact  point  on  the  ground  plane,  $  *  o  is  defined  as  the  orientation  of 
maximum  velocity  on  the  ground  plane.  Donaldson  and  Snedeker,  Reference  13, 
have  shown  that  for  oblique  jet  impingement,  the  maximum  pressure  point  on  the 
ground  plane  is  displaced  away  from  the  jet  centerline  in  the  direction 
<(>  -  180°.  The  shifted  stagnation  point  represents  the  effective  origin  of  the 
wall-jet  flow  that  is  described  by  the  coordinates  r  and  The  magnitude  of 
this  displacement,  denoted  by  AR^,  depends  on  both  jet  inclination  angle  and 
distance  from  the  jet  exit  to  the  ground.  In  the  present  analysis,  ARj./6^  as 

a  function  of  is  computed  using  the  curve  of  Donaldson  and  Snedeker  which 

is  reproduced  in  Figure  18(c). 


With  the  origin  of  the  wall  jet  established,  the  wall-jet  velocity 

profile  is  defined  using  the  methodology  presented  in  Reference  1.  The 

maximum  wall-jet  velocity  at  the  edge  of  the  impingement  region,  U  , 

max 

is  related  to  V.  by 
JCI 
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U  =  V .  0.55  +  0.170 

max].  jCj. 


x  1.0  -  0.4 


/90°  -  °1lV 

0.170  - ^ 

\  15°  / J 


V.  is  determined  from  the  curve  of  Figure  18(b).  The  behavior  of  Umax  with 
JCI 

r  follows  the  characteristic  radial  decay, 


U  (r)  -  —  U 
max  r  raax^ 


(A-7) 


where  Rj  is  evaluated  from  the  correlation 


Rt  «  4.8  fij. 


Sj  follows  from  the  curve  of  Figure  18(b).  Knowing  13max  at  the  required 
points  on  the  outflow  boundaries,  a  polynomial  curve  fit  to  the  characteristic 
wall-jet  velocity  profile  shown  in  Figure  18(d)  is  used  to  compute  the 
velocity  profile  as  a  function  of  N/Nj.  N  denotes  the  normal  distance  above 
the  ground  plane,  and  N5  denotes  the  value  at  which  U  “0.5  Umax*  N5,  a 
function  of  r,  is  computed  from 


N5  (r)  -  N5  +  0.07  (r  -  Rj), 


with  evaluated  from  a  mass  balance  given  in  Reference  1.  Equations  (A-6) 

5 1 

through  (A-9)  are  used  to  compute  the  wall-jet  velocity  profiles  resulting 
from  each  of  the  primary  jets,  and  the  projections  of  U  at  the  required  points 
on  the  outflow  planes  are  computed. 

The  final  step  is  to  compute  the  velocity  profiles  through  the  fountain 
on  the  computational  domain  boundaries  y  -  h  and  x  *  w.  This  step  requires 
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calculation  of  the  intersection  of  the  fountain  with  these  planes,  a 
computation  which  begins  with  the  determination  of  the  fountain  stagnation¬ 
line  displacement  on  the  ground  plane.  The  distance  from  the  origin  of  wall 
jet  1  to  the  point  where  the  opposing  wall  jets  collide  is  given  by 


r  =  AR  . 
l,s  I, l 


2  sin  a 


+  z 


JI 


(  A- 1 0 ) 


Similarly,  the  distance  from  the  origin  of  wall  jet  2  to  the  collision  line  is 
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2,s 
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2  sin  a . T 

Jl 


(A-l 1 ) 


In  the  previous  equations,  zq'  represents  the  displacement  of  the  fountain 
stagnation  line  from  the  midpoint  of  the  distance  between  the  intercepts  of 
the  geometric  nozzle  centerlines  and  the  ground  plane.  Using  Equation  (A— 7 ) 
and  the  condition  that  the  wall-jet  velocities  are  equal  at  the  collision 
line,  it  follows  that 
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max 
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max 


(A -12) 


The  combination  of  Equations  (A-10)  through  (A-12)  gives 
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(A-14) 
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1,1  2  sin  a 
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and 


f2(ajl) 


2  sin  a 


-  AR 


jl 


1,2 


(A-15) 


Consider  the  velocity  ratio  U  /U  in  Equation  (A-13).  Wall  jet  1 

J  max  „  max  , 

L  f  Z  1  9 

approaches  the  stagnation  line  with  (j>  =  0,  and  wall  jet  2  approaches  the 
stagnation  line  with  4>  =  180°.  Therefore,  from  Equation  (A-6), 


(A-l 6) 


Equation  (A-14)  in  conjunction  with  Equation  (A-13)  provides  the  displacement 
of  the  stagnation  line. 

The  next  step  is  to  compute  the  inclination  of  the  fountain.  Using 
conservation  of  monemtum  principles  and  the  approximation  that  the  upwash  can 
be  treated  as  a  radial  planar  flow,  it  can  be  shown  in  an  analysis  analogous 
to  that  given  in  Reference  10  that  the  upwash  inclination  0  is  related  to  the 
ratio  of  the  wall-jet  thicknesses  in  the  following  manner: 


0  »  sin 


(A-l 7 ) 


where 
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With  the  geometry  of  the  fountain  established,  the  velocity  profiles  can 
be  computed  on  planes  passing  through  the  fountain  which  are  the  boundaries  of 
the  computational  domain.  Following  the  analysis  of  Jenkins  and  Hill, 
Reference  10,  the  maximum  upwash  velocity  Wfflax  at  a  specified  value  of  y  for 
x  =  0  is  given  by 


(A-19) 


where  Pg  is  the  surface  pressure  at  the  stagnation  line,  P^  is  the  ambient 
pressure,  and  j  is  the  dynamic  pressure  at  the  jet  exit.  The  pressure 
coefficient  term  is  evaluated  from  an  empirical  correlation.  The  maximum 
local  velocity  in  the  upwash  is  assumed  to  have  the  following  variation  with 
radial  distance  from  the  virtual  origin: 
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W  (r  )  =  —  W 
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where 
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Therefore, 
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w  (r  ) 
max  u 
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max 


AR  .  +  -z - : -  +  z  '  +  y' 

1,1  2  sin  a . T  o  ' 

_ J1 
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Finally,  the  velocity  variation  normal  to  the  upwash  plane,  w(ru>  £)»  is 
computed,  where  £  denotes  a  coordinate  perpendicular  to  this  plane  passing 
through  ru.  The  assumption  is  made  that  the  profile  has  the  Gaussian  form 
given  in  Reference  15, 


"  <'„•  S)  ,  2, 

H  (r  )  '  “P  [-0.693(E/{0-25)  ]. 

max  u 


(A-24) 


In  this  equation  F  represents  the  distance  from  the  center  of  the  upwash 

U  •  Z  Z) 

to  the  half-velocity  point  of  the  profile  and  is  computed  from  the  following 
relationship : 
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Substitution  of  Equation  (A-24)  into  Equation  (A-23)  yields 


ii  i  n  (ARt  ,  +7 — 7^ -  +z  '  +  y'l  exp  [-0.693  (£  /F  )  "  1 

W  (ry;  K)  y  1,1  2  sin  o  J  1  0.25 
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I  1,1  2  sin  opjj  o  J  J 
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(A—  26) 


The  required  components  of  the  velocity  profiles  normal  to  the  outflow  planes 
are  obtained  using  simple  geometric  relationships.  The  velocity  variation 
over  the  remainder  of  the  top  surface  of  the  computational  domain  is  computed 
using  an  integral  mass  balance  in  which  the  difference  between  the  outflow  and 
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inflow  of  mass  is  taken  to  be  the  mass  entrained.  The  latter  is  distributed 
uniformly  outside  the  primary  jets  and  fountain  and  is  used  to  evaluate  the 
entrainment  velocity  distribution. 

The  specification  of  the  boundary  conditions  is  simplified  for  the  case 
of  two  jets  with  normal  impingement.  In  this  configuration  (illustrated  in 
Figure  2-(a)),  two  additional  symmetry  planes  comprise  the  boundaries  of  the 
computational  domain,  and  the  wall-jet  velocity  distribution  need  be 
calculated  on  only  a  single  outflow  plane. 

It  should  be  emphasized  that  the  numerical  solution  scheme  for  the  time- 
averaged  Navier-Stokes  equations  described  in  Section  3.0  is  not  limited  to 
the  boundary  conditions  given  in  this  appendix.  As  improvements  are  made  in 
the  engineering  approximations  used  to  define  the  boundary  values  of  the 
normal  velocity  components,  the  corresponding  changes  in  boundary  conditions 
need  be  made  in  only  one  routine  of  the  computer  code. 
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APPENDIX  B:  CENTRAL-DIFFERENCE  FORMULAS  USED  IN  THE  FINITE-DIFFERENCE 
SOLUTION  OF  THE  GOVERNING  EQUATIONS 


In  the  numerical  solution  scheme  described  in  Section  3.0,  central- 
difference  approximations  to  the  following  derivatives  of  (j>  (an  arbitrary  flow 
variable)  are  required: 


and 

3x3z 


The  approximate  forms  of  these  derivatives  are  derived  using  a  Taylor- 
series  expansion  and  retaining  a  sufficient  number  of  terms  to  ensure  that 
each  approximating  formula  is  accurate  to  0(h),  where  h  is  the  finite- 
difference  node  spacing  in  all  three  coordinate  directions.  Different 
formulas  are  required  for  the  interior  region,  the  boundary  planes,  and  the 
corners  of  the  computational  domain.  Representative  formulas  for  derivatives 
with  respect  to  x  are  given  below  for  the  interior,  the  boundary  plane  i  =  1, 
and  the  lower  corner  i  *  1,  j  «*  1  of  the  nodal  network.  Analogous  formulas 
apply  to  the  boundary  plane  i  «  I  and  the  remaining  corners  and  to  the 
derivatives  with  respect  to  y  and  z.  The  notation  used  in  the  following 
finite-difference  approximations  is  defined  in  Figure  3. 
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3"i> 
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(corner  point,  i  *  1  and  j  «  1) 


The  derivative  formulas  are  stored  on  a  library  tape  which  is  accessed  by 
the  computer  codes  used  to  solve  the  finite-difference  equations. 
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